---
title: "Ukraine & Ethiopia"
output: pdf_document
date: "2023-02-05"
---


```{r include=FALSE}
rm(list=ls())
```

```{r eval=FALSE, include=FALSE}
install.packages("ggplot2")
install.packages("ggmap")
install.packages("rgdal")
install.packages("rgeos")
install.packages("maptools")
install.packages("dplyr")
install.packages("tidyr")
install.packages("knitr")
install.packages("spData")
install.packages("Rtools")
install.packages("tidyverse")
install.packages("lmtest")
install.packages("plm")
install.packages("aod")
install.packages("cartogram")
install.packages("broom")
install.packages("viridis")
install.packages("patchwork")
install.packages("geojsonio")
install.packages("sjmisc")
install.packages("interactions")
install.packages("sjPlot")
install.packages(("gtable"))
install.packages('grid')  
install.packages("parsedate")
install.packages("splancs")
install.packages("st")
install.packages('osmdata')
install.packages("ggpubr")
install.packages("spatstat")
install.packages("doParallel")
install.packages("meltt")
install.packages("geosphere")
install.packages("nngeo")
install.packages("ggspatial")
library(ggspatial)

```

```{r include=FALSE}

knitr::opts_chunk$set(echo = FALSE)
require(ggplot2)
require(rgdal)
require(maptools)
require(rgeos)
require(dplyr)
require(tidyr)
library(leaflet)
library(raster)
library(dplyr)
library(spData)
library(mapview)
library(shiny)
library(tidyverse)
library(devtools)
library(cartogram)
library(broom)
library(viridis)
library(geojsonio)
library(sjmisc)
library(sandwich)
library(texreg)
library(interactions)
library(RColorBrewer)
library(sp)
library(readxl)
library(ggthemes)
library(splancs) 
library(parsedate)
library(st)
library(gridExtra)
library(ggmap)
library(osmdata)
library(patchwork)
library(spatstat)
library(rgeos)
library(sf)
library(doParallel)
library(foreach)
library(meltt)
library(geosphere)
library(spatstat)
library(nngeo)
library(spdep)
library(ggmap)
library(ggsn)

```

    
```{r include=FALSE}
 
# UCDP Candidate Events Dataset (UCDP Candidate) version 22.01.22.12
geo<-"https://ucdp.uu.se/downloads/candidateged/GEDEvent_v22_01_22_12.csv"
download.file(geo, basename(geo), mode="wb")
geomont<-read.csv("GEDEvent_v22_01_22_12.csv")
ukr<-subset(geomont,geomont$country=="Ukraine")

```

```{r include=FALSE}
#ACLED for Ukraine
setwd("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr")
acled<-read.csv("acled.csv")

#Ukraine
ukraine<-subset(acled,acled$country=="Ukraine" & year=="2022")

#East
ukrainewar<-subset(acled,acled$admin1=="Luhansk"|
                     acled$admin1=="Kharkiv"|
                     acled$admin1=="Donetsk"|
                     acled$admin1=="Zaporizhia"|
                     acled$admin1=="Dnipropetrovsk")

```

```{r}
#Adjust Dates and limit to explosions
ukrainewar$event_date<-parse_date(ukrainewar$event_date)
ukrainewar$event_date<-as.Date(ukrainewar$event_date)

#Date range
ukrainewar<-subset(ukrainewar, ukrainewar$event_date >="2021-11-26" &
                     ukrainewar$event_date <= "2022-5-25")

ukrainepre<-subset(ukrainewar,ukrainewar$event_date <= "2022-02-23")

ukrainepost<-subset(ukrainewar,ukrainewar$event_date >= "2022-02-24")

```


```{r include=FALSE}
#Ukraine Map

#admin level 1
ukrmap1<-"https://stacks.stanford.edu/file/druid:gg870xt4706/data.zip"
download.file(ukrmap1,basename(ukrmap1),mode="wb")
unzip("data.zip")
ukraine_poly<-readOGR("UKR_adm1.shp")
ukraine<-read.csv("UKR_adm1.csv")

write.csv

#admin level 2
ukrmap2<-"https://stacks.stanford.edu/file/druid:pp624tm0074/data.zip"
download.file(ukrmap2,basename(ukrmap2),mode="wb")
unzip("data.zip")
ukra2<-readOGR("UKR_adm2.shp")
```


```{r include=FALSE}
#Limit Geographically

unique(ukraine_poly$NAME_1)

ukrainepl<-subset(ukraine_poly,
                  ukraine_poly$NAME_1=="Luhans'k"|
                    ukraine_poly$NAME_1=="Donets'k"|
                    ukraine_poly$NAME_1=="Kharkiv"|
                    ukraine_poly$NAME_1=="Dnipropetrovs'k"|
                    ukraine_poly$NAME_1=="Zaporizhzhya")

plot(ukrainepl)

ukraine_df<-fortify(ukrainepl)

``` 


```{r}
#Patterns of War
ukra2_df<-fortify(ukra2)
ukraine_df$ind<-1

p2<-ggplot()+
  geom_path(data=ukra2_df,aes(long,lat,group=group))+
  geom_polygon(data=ukraine_df,aes(fill=ind,long,lat,group=group),
               color="black",show.legend=FALSE)+ 
  geom_point(data=ukrainepre, aes(x=longitude, y=latitude), 
             color="black", size=1)

p2<-p2+theme_map()

p3<-ggplot()+
  geom_path(data=ukra2_df,aes(long,lat,group=group))+
  geom_polygon(data=ukraine_df,aes(fill=ind,long,lat,group=group),
               color="black",show.legend=FALSE)+
  geom_point(data=ukrainepost, aes(x=longitude, y=latitude),
             color="black", size=1)

p3<-p3+theme_map()

grid.arrange(p2, p3, nrow = 1) 

```


 
```{r include=FALSE}
#Firedata for Ukraine

#NOAA-20
ukrfire1<-read.csv ("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/fire_nrt_J1V-C2_358261.csv")

#Suomi Not usde
ukrfire2<-read.csv ("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/fire_archive_SV-C2_358262.csv")

 
 
```

```{r include=FALSE}
#Reformat Date
ukrfire1$acq_date<-as.Date(ukrfire1$acq_date)
class(ukrfire1$acq_date)


#90 days before and after 
ukrfire1<-subset(ukrfire1,ukrfire1$acq_date >= "2021-11-26" 
                 & ukrfire1$acq_date <= "2022-05-25")
ukrfirepre<-subset(ukrfire1,ukrfire1$acq_date <= "2022-02-23")
ukrfirepost<-subset(ukrfire1,ukrfire1$acq_date >= "2022-02-24")

```

```{r include=FALSE}
#Plot
f1<-ggplot()+
  geom_path(data=ukra2_df,aes(long,lat,group=group))+
  geom_polygon(data=ukraine_df,aes(fill=ind,long,lat,group=group),
               color="black",show.legend=FALSE)+
  geom_point(data=ukrfirepre, aes(x=longitude, y=latitude), color="black", size=1)

f1<-f1+theme_map()
f1

f2<-ggplot()+geom_path(data=ukra2_df,aes(long,lat,group=group))+
  geom_polygon(data=ukraine_df,aes(fill=ind,long,lat,group=group),
               color="black",show.legend=FALSE)+
  geom_point(data=ukrfirepost, aes(x=longitude, y=latitude),
             color="black", size=1)

f2<-f2+theme_map()
f2


```

```{r}
crs <- st_crs("+proj=longlat +datum=WGS84 +no_defs")
violence <- ukrainepost
```


```{r eval=FALSE, include=FALSE}
#Prepare files for ARCGIS.

fires <- ukrfirepost
fires <- st_as_sf(fires, coords = c("longitude", "latitude"))
fires <- as(fires, "sf")
fires <- st_set_crs(fires, crs)
fires <- st_transform(fires, crs)

st_write(fires, "C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/ARCGIS/fires.shp")


firespre <- ukrfirepre
firespre <- st_as_sf(firespre, coords = c("longitude", "latitude"))
firespre <- as(firespre, "sf")
firespre <- st_set_crs(firespre, crs)
firespre <- st_transform(firespre, crs)

st_write(firespre, "C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/ARCGIS/firespre.shp")

violence <- st_as_sf(violence, coords = c("longitude", "latitude"))
violence <- as(violence, "sf")
violence <- st_set_crs(violence, crs)
violence <- st_transform(violence, crs)

st_write(violence, "C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/ARCGIS/violence.shp")

#As csv
write_csv(ukrfirepost, "C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/ARCGIS/NOAAfires.csv")

write_csv(ukrainepost, "C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/ARCGIS/ukrainepost.csv")


```


```{r}
#admin
ukrasf<-read_sf("UKR_adm2.shp")
ukrasf <- subset(ukrasf, NAME_1 %in% 
                   c("Luhans'k", "Donets'k", 
                     "Kharkiv", "Dnipropetrovs'k", "Zaporizhzhya"))

ukrasf <- st_set_crs(ukrasf, crs)

```

```{r eval=FALSE, include=FALSE}
st_write(ukrasf, "C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/admin2.shp")
```


```{r}
#Fires as data frames
eastfirespre<-read.csv("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/Eastfirespre.csv")

eastfirespost<-read.csv("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/Eastfirespost.csv")

```


```{r}
#Fires and ACLED density
proj4string(ukrainepl) = 
  CRS("+proj=utm +zone=46 +datum=WGS84 +units=m +no_defs")


coordinates(eastfirespre) <- ~Long+Lat
ukrOwin <- as.owin(ukrainepl)
proj4string(eastfirespre) = 
  CRS("+proj=utm +zone=46 +datum=WGS84 +units=m +no_defs")
pts1<-coordinates(eastfirespre)


coordinates(eastfirespost) <- ~Long+Lat
ukrOwin <- as.owin(ukrainepl)
proj4string(eastfirespost) = 
  CRS("+proj=utm +zone=46 +datum=WGS84 +units=m +no_defs")
pts2<-coordinates(eastfirespost)

#Conflict
coordinates(ukrainepre) <- ~longitude+latitude
proj4string(ukrainepre) = 
  CRS("+proj=utm +zone=46 +datum=WGS84 +units=m +no_defs")
pts3<-coordinates(ukrainepre)

  
coordinates(ukrainepost) <- ~longitude+latitude
proj4string(ukrainepost) = 
  CRS("+proj=utm +zone=46 +datum=WGS84 +units=m +no_defs")
pts4<-coordinates(ukrainepost)

```


```{r echo=TRUE}
#Point Pattern Object
part1 <- ppp(pts1[,1], pts1[,2], window=ukrOwin)
part2 <- ppp(pts2[,1], pts2[,2], window=ukrOwin)
part3 <- ppp(pts3[,1], pts3[,2], window=ukrOwin)
part4 <- ppp(pts4[,1], pts4[,2], window=ukrOwin)

ds1 <- density(part1,bw=1,adjust=0.4)
ds2 <- density(part2,bw=1,adjust=0.2)
ds3 <- density(part3,bw=1,adjust=0.2)
ds4 <- density(part4,bw=1,adjust=0.2)

#COlours
col1 <-  colorRampPalette(c("black", "green", "red","orange"))(100)

#To highlight DS1
ds1Broad <- density(part1,bw=1,adjust=0.4)

plot(ds1Broad, main = "Fire Density Pre-Invasion")
```

```{r echo=TRUE}
 
layout(matrix(c(1, 3, 2, 4), nrow = 2, byrow = TRUE))

par(mar = c(0.5, 0.25, 0.75, 2) + 0.1, oma = c(0, 0, 2, 0))

plot(ds1, main = "Fire Density Pre-Invasion",col = col1)
plot(ds2, main = "Fire Density Post-Invasion",col = col1)
plot(ds3, main = "ACLED Density Pre-Invasion",col = col1)
plot(ds4, main = "ACLED Density Post-Invasion",col = col1)



```



```{r include=FALSE}
 
# Firedata in Relevant Region
firestable <- read.csv("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/NOAAFirestable.csv") 
firestableSUOMI <- read.csv("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/Firestable.csv") 
violence2 <- ukrainepost
violence2 <- data.frame(violence2)

firestable$type_tax <- "fire"
violence2$type_tax <- "conflict"

firestable$acq_date <- as.Date(firestable$acq_date)
violence2$event_date <- as.Date(violence2$event_date)

firestable <- subset(firestable, select = -c(Long, Lat,Join_Count,TARGET_FID))

 
 
```


  
```{r}
#ACLED intersected with Admin 2 for regional matching
violadmin <- read.csv("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/intadmin2violence.csv")
violadmin$event_dt <- as.Date(violadmin$evnt_dt)
firestable$acq_date <- as.Date(firestable$acq_date)
names(violadmin)[names(violadmin) == "Long"] <- "longitude"
names(violadmin)[names(violadmin) == "Lat"] <- "latitude"
names(violadmin)[names(violadmin) == "Name_2222"] <- "NAME_2"
names(firestable)[names(firestable) == "acq_date"] <- "date"
names(violadmin)[names(violadmin) == "evnt_dt"] <- "date"

violadmin$date <- as.Date(violadmin$date)
firestable$date <- as.Date(firestable$date)

```
 
```{r}
# Create an index for faster subsetting
violadmin_index <- with(violadmin, paste(NAME_2, date))

# Initialize an empty Match_Status vector
match_status <- character(nrow(firestable))

# Iterate over each row in firestable
for (i in 1:nrow(firestable)) {
  # Get the current event NAME_2 and date from firestable
  event_NAME_2 <- firestable$NAME_2[i]
  event_date <- as.Date(firestable$date[i])
  
  # Check if the combination of NAME_2 and date exists in violadmin
  index <- paste(event_NAME_2, event_date)
  if (index %in% violadmin_index) {
    match_status[i] <- "Match"
  } else {
    match_status[i] <- "Non-Match"
  }
}

# Add the Match_Status column to firestable
matchedregion <- data.frame(firestable, Match_Status = match_status, stringsAsFactors = FALSE)

```

 

  

```{r}
# Set the maximum distance in meters and time difference in seconds
max_distance <- 5000
max_time_diff <- 24 * 60 * 60  # 1 day in seconds

# Empty dataframe to store matched rows
matcheddist <- data.frame()

# Iterate over each row in firestable
for (i in 1:nrow(firestable)) {
  # Get the current event coordinates and timestamp from firestable
  event_lat <- firestable$latitude[i]
  event_lon <- firestable$longitude[i]
  event_time <- as.POSIXct(firestable$date[i])
  
  # Filter violence2 dataframe based on distance and time difference criteria
  filtered_rows <- with(violence2, distHaversine(matrix(c(event_lon, event_lat), ncol = 2, byrow = TRUE), cbind(longitude, latitude)) <= max_distance &
                                   abs(as.numeric(difftime(event_time, as.POSIXct(event_date), units = "secs"))) <= max_time_diff)
  
  # Check for match 
  if (any(filtered_rows)) {
    matcheddist <- rbind(matcheddist, cbind(firestable[i, ], Match_Status = "Match"))
  } else {
    matcheddist <- rbind(matcheddist, cbind(firestable[i, ], Match_Status = "Non-Match"))
  }
}
```

```{r eval=FALSE, include=FALSE}

range(firestable$longitude)
range(firestable$latitude)

range(violadmin$longitude)
range(violadmin$latitude)
```



```{r}
#The difference of K-functions with Monte Carlo envelopes.
 
window1 <- owin(xrange = c(32, 41), yrange = c(46, 51))
firestable_pp <- ppp(firestable$longitude, firestable$latitude, window = window1)

window2 <- owin(xrange = c(32 , 41), yrange = c(46 , 51))
violadmin_pp <- ppp(violadmin$longitude, violadmin$latitude, window = window2)


K_firestable <- Kest(firestable_pp,correction="Ripley")
K_violadmin <- Kest(violadmin_pp,correction="Ripley")

```

```{r}
# With Montecarlo Simulation

# Generate Monte Carlo envelopes for the difference
envelope_diff1 <- envelope(firestable_pp, Kest, correction = "Ripley", verbose = F)

envelope_diff2 <- envelope(violadmin_pp, Kest, correction = "Ripley", verbose = F)

#Figure 2 in the Appendix
plot(K_firestable, main = "", xlab = "r", ylab = "K(r)", col = "red", xlim = c(0, max(K_firestable$r)))
plot(envelope_diff1,main = "")

#Figure 3 in the Appendix
plot(K_violadmin, main = "",  xlab = "r", ylab = "K(r)", col = "blue", xlim = c(0, max(K_violadmin$r)))
plot(envelope_diff2,main = "")

```

```{r}
#L Function
envelope_diff3 <- envelope(firestable_pp, Lest, correction = "Ripley", verbose = F)
plot(envelope_diff3, . - r ~ r,main = "")

envelope_diff4 <- envelope(violadmin_pp, Lest, correction = "Ripley", verbose = F)
plot(envelope_diff4, . - r ~ r,main = "")
```

```{r}
#inhomogenous test
#L-hat estimates for inhomogeneous Poisson point processes with envelopes adjusted for density
plot(envelope(firestable_pp, Linhom, sigma = 0.75, correction = "Ripley", verbose = F), 
    . - r ~ r, main = "")
plot(envelope(violadmin_pp, Linhom, sigma = 0.75, correction = "Ripley", verbose = F), 
    . - r ~ r, main = "")


```

```{r}
#Cross-K
K_diff <- Kest(firestable_pp, violadmin_pp)

# Generate Monte Carlo envelopes for the difference
envelope_diff <- envelope(firestable_pp, Kest, nsim = 99, sim = "independent", nsimE = 99, correction = "border", verbose = FALSE)

plot(K_firestable, main = " ", xlab = "r", ylab = "K_diff(r)")
plot(envelope_diff, type = "envelope", col = "lightgray", lwd = 2, add = TRUE)

```


```{r}
#Pairwise correlation function
K <- Kest(firestable_pp, violadmin_pp, correction = "iso")
g <- pcf(firestable_pp, violadmin_pp, correction = "iso")
plot(anylist(K, g), main = "", main.panel = "")
 
# Convert KM to meters
dist_meters <- K$dist * 1000

plot(g, main = "", main.panel = "", xlab = "Distance (meters)")
axis(side = 1, at = dist_meters, labels = dist_meters)

```
 
 


```{r}
 
# Study region encompassing both datasets
bb <- boundingbox(c(firestable_pp, violadmin_pp))
W <- as.owin(bb)

# Combine 
combined_pp <- superimpose(firestable_pp, violadmin_pp, W=NULL)

# FitStrauss process model
model <- ppm(combined_pp, ~1, Strauss(r = 0.1))

summary(model)


```

```{r}
plot(firestable_pp, pch=16)
plot(violadmin_pp, pch=16)

```
 
```{r}
#Buffered Violence Points
viol0km <- st_read("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/violence0/Violence.shp")
viol1km <- st_read("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/violence1km/BufferedFeatures.shp")
viol3km <- st_read("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/violence3km/BufferedFeatures.shp")
viol5km <- st_read("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/violence5km/BufferedFeatures.shp")
viol10km <- st_read("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/violence10km/BufferedFeatures.shp")
viol20km <- st_read("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/violence20km/BufferedFeatures.shp")
fireshp <- st_read("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/fires_shp/OverlayedFeatures.shp")

```

```{r}
#adjust date format
fireshp$date <- as.Date(fireshp$acq_date, format = "%Y-%m-%d")
viol0km$date <- as.Date(viol0km$evnt_dt, format = "%Y-%m-%d")
viol1km$date <- as.Date(viol1km$evnt_dt, format = "%Y-%m-%d")
viol3km$date <- as.Date(viol3km$evnt_dt, format = "%Y-%m-%d")
viol5km$date <- as.Date(viol5km$evnt_dt, format = "%Y-%m-%d")
viol10km$date <- as.Date(viol10km$evnt_dt, format = "%Y-%m-%d")
viol20km$date <- as.Date(viol20km$evnt_dt, format = "%Y-%m-%d")

```

```{r warning=FALSE, include=FALSE}
#Intersects on the same day

#0KM
fireshp$intersects_viol0km <- st_intersects(fireshp, viol0km) %>%
  apply(1, function(x) any(viol0km$date == fireshp$date[x]))

num_intersect0 <- sum(fireshp$intersects_viol0km)

#1KM
fireshp$intersects_viol1km <- st_intersects(fireshp, viol1km) %>%
  apply(1, function(x) any(viol1km$date == fireshp$date[x]))

num_intersect1 <- sum(fireshp$intersects_viol1km)

#3km
fireshp$intersects_viol3km <- st_intersects(fireshp, viol3km) %>%
  apply(1, function(x) any(viol3km$date == fireshp$date[x]))

num_intersect2 <- sum(fireshp$intersects_viol3km)

#5km
fireshp$intersects_viol5km <- st_intersects(fireshp, viol5km) %>%
  apply(1, function(x) any(viol5km$date == fireshp$date[x]))

num_intersect3 <- sum(fireshp$intersects_viol5km)

#10km
fireshp$intersects_viol10km <- st_intersects(fireshp, viol10km) %>%
  apply(1, function(x) any(viol10km$date == fireshp$date[x]))

num_intersect4 <- sum(fireshp$intersects_viol10km)

#20KM
fireshp$intersects_viol20km <- st_intersects(fireshp, viol20km) %>%
  apply(1, function(x) any(viol20km$date == fireshp$date[x]))

num_intersect5 <- sum(fireshp$intersects_viol20km)

```

```{r}
#Percentages for Figure 4 
probability0 <- num_intersect0 / nrow(fireshp)
probability0
probability1 <- num_intersect1 / nrow(fireshp)
probability1
probability2 <- num_intersect2 / nrow(fireshp)
probability2
probability3 <- num_intersect3 / nrow(fireshp)
probability3
probability4 <- num_intersect4 / nrow(fireshp)
probability4
probability5 <- num_intersect5 / nrow(fireshp)
probability5

 
```

```{r}
#Figure 4 

#Vector of X-axis values
x <- c(1, 3, 5, 10, 20)

#Vector of Y-axis values
y <- c(0.082, 0.382, 0.606, 0.827, 0.918)

# Plot
plot(x, y, type = "b", pch = 16, ylim = c(0, 1), xlab = "Buffer Radius in km", ylab = "Proportion of Intersections", xaxt = "n")
grid()
title(main = "")
text(x, y, labels = y, pos = ifelse(x == 1, 4, 3))
axis(1, at = x, labels = x)
 
 
```
 
 
 
```{r}
# Step 1: Filter events occurring on the same day
same_day_data <- inner_join(violadmin, firestable, by = "date")

# Step 2: Calculate the distance between observations
same_day_data <- same_day_data %>%
  mutate(distance = distVincentySphere(cbind(longitude.x, latitude.x),
                                       cbind(longitude.y, latitude.y)))

# Step 3: Group by the observation in the violadmin dataset
# and select the minimum distance
nearest_neighbors <- same_day_data %>%
  group_by(longitude.x, latitude.x) %>%
  summarise(min_distance = min(distance))

# Step 4: View the results
print(nearest_neighbors)

```

```{r}
# Step 1: Filter events occurring on the same day
same_day_data <- inner_join(violadmin, firestable, by = "date")

# Step 2: Calculate the distance between observations
same_day_data <- same_day_data %>%
  mutate(distance = distVincentySphere(cbind(longitude.x, latitude.x),
                                       cbind(longitude.y, latitude.y)))

# Step 3: Group by the observation in the violadmin dataset
# and select the k nearest neighbors
k_nearest <- 10  # Set the desired maximum number of nearest neighbors

average_distances <- data.frame(k = 1:k_nearest, avg_distance = numeric(k_nearest))

for (k in 1:k_nearest) {
  # Find k nearest neighbors for each observation
  nearest_neighbors <- same_day_data %>%
    group_by(longitude.x, latitude.x) %>%
    top_n(k, -distance) %>%
    ungroup()

  # Calculate the average distance for k nearest neighbors
  avg_distance <- mean(nearest_neighbors$distance)

  # Store the average distance in the data frame
  average_distances[k, "avg_distance"] <- avg_distance
}

# Step 4: View the results
print(average_distances)

```

```{r}
# Figure 5
ggplot(average_distances, aes(x = k, y = avg_distance)) +
  geom_line() +
  labs(x = "Number of Nearest Neighbors", y = "Average Distance (meters)") +
  ggtitle("") +
  scale_x_continuous(breaks = seq(1, 20, by = 1)) +
  theme(plot.margin = margin(30, 30, 30, 30), 
        axis.title.x = element_text(margin = margin(t = 10)),
        axis.title.y = element_text(margin = margin(r = 10)))


```

```{r}
#Appendix on undetected events in Ukraine

#Populated areas, polygons: https://data.humdata.org/dataset/hotosm_ukr_populated_places

#I use ARCGIS online due to the significant computing power required.  

#populated areas in the ROIs
pop_rois<- st_read("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/populated_polygons/OverlayedFeatures.shp")
base_plot<-ggplot() + geom_sf(data = pop_rois) +   theme_void()
 
```


```{r}

# Perform a spatial join between "fireshp" and "pop_rois" based on spatial intersection
fireshp_intersect <- st_intersection(fireshp, pop_rois)

# Perform a spatial join between "violadmin" and "pop_rois" based on spatial intersection
viol0km_intersect <- st_intersection(viol0km, pop_rois)

base_plot<-ggplot() + geom_sf(data = pop_rois) +   theme_void()
base_plot

plot <- ggplot() +
  geom_sf(data = pop_rois) +
  geom_sf(data = fireshp_intersect, aes(color = "Fires")) +
  geom_sf(data = viol0km_intersect, aes(color = "Conflict Events")) +
  scale_color_manual(
    values = c("blue", "red"),
    labels = c("Conflict events", "Fires"),
    guide = guide_legend(title = "Event Types")
  ) +
  theme_void()

# Display the plot
print(plot)

```

```{r}
fireshp_intersect <- st_intersection(fireshp, pop_rois)

viol0km_intersect <- st_intersection(viol0km, pop_rois)

base_plot<-ggplot() + geom_sf(data = pop_rois) +   theme_void()
base_plot

plot <- ggplot() +
  geom_sf(data = pop_rois) +
  geom_sf(data = fireshp_intersect, aes(color = "Fires")) +
  geom_sf(data = viol0km_intersect, aes(color = "Conflict Events")) +
  scale_color_manual(
    values = c("red", "blue"),
    labels = c("Fires", "Conflict Events"),
    guide = guide_legend(title = "Event Types")
  ) +
  theme_void()

```
 
```{r eval=FALSE, include=FALSE}
#Prepare for ArCGIS
st_write(fireshp_intersect, "C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/ARCGIS/fireshp_intersect.shp")
st_write(viol0km_intersect, "C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/ARCGIS/viol0km_intersect.shp")

```

```{r}
#Fire detections not in shared polygon
firenoviol<- st_read("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/firenonmatch/fireshp_intersect_selection__2_.shp")
 
```
```{r}

plot <- ggplot() +
  geom_sf(data = pop_rois) +
  geom_sf(data = fireshp_intersect, aes(color = "Fires")) +
  geom_sf(data = viol0km_intersect, aes(color = "Conflict Events")) +
  geom_sf(data = firenoviol, aes(color = "Non-Match Fires")) +
  scale_color_manual(
    values = c("blue", "red", "green"),
    labels = c("Conflict Events", "Matched Fires", "UnmatchedS Fires"),
    guide = guide_legend(title = "Event Types")
  ) +
  theme_void()

plot


```
```{r eval=FALSE, include=FALSE}
#Random Sample
set.seed(123)
sample_size <- 100
random_sample <- firenoviol[sample(nrow(firenoviol), sample_size, replace = FALSE), ]

 
```

```{r eval=FALSE, include=FALSE}
#Sample Dataset as CSV
rs <- subset(random_sample, select = c("acq_dat", "geometry","VARNAME_2"))
rs_sf <- st_as_sf(rs, coords = c("longitude", "latitude"), crs = 4326)
rs_sf <- st_transform(rs_sf, crs = 4326)
rs$latitude <- st_coordinates(rs_sf)[, "Y"]
rs$longitude <- st_coordinates(rs_sf)[, "X"]
rs$geometry <- NULL
rs$ID <- sprintf("FIR-%03d", seq_len(nrow(rs)))
rs$prob<-""
rs$URL <- ""
rs$note <- ""
colnames(rs) <- c("acq_dat", "NAME", "latitude", "longitude",   "ID", "prob","URL", "note")

write.csv(rs, "C:/Users/Mikae/OneDrive/Dokumenter/Articles/Firedata Ukr/H3.csv")

```



























```{r}
#Appendix Table 1 

#matched data types
matcheddist2 <- subset(matcheddist, matcheddist$Match_Status=="Match")

type_table1 <- table(matcheddist2$violence2)
typeprop1 <- prop.table(type_table1) * 100

rt1 <- data.frame(Event_Type = names(type_table1),
                           Percentage = typeprop1)

print(rt1) 
 
```

```{r}
#Figure 1 in the Appendix

# Counts and proportions for fire events
fires_counts <- ukrfire1 %>%
  group_by(acq_date) %>%
  summarise(count = n()) %>%
  mutate(proportion = count / sum(count))

# Counts and proportions for violence events
violence_counts <- ukrainewar %>%
  group_by(event_date) %>%
  summarise(count = n()) %>%
  mutate(proportion = count / sum(count))

fires_counts$date <- as.Date(fires_counts$acq_date)
violence_counts$date <- as.Date(violence_counts$event_date)

ggplot() +
  geom_line(data = fires_counts, aes(x = date, y = proportion, color = "Fires")) +
  geom_line(data = violence_counts, aes(x = date, y = proportion, color = "Violence")) +
  labs(x = "", y = "Proportion of Observations") +
  scale_x_date(date_breaks = "1 month", date_labels = "%b %Y") +
  scale_color_manual(values = c("Fires" = "red", "Violence" = "blue"))

```

```{r}
#Random Sample
firesampl<-data.frame(fireshp)
set.seed(123)
```
  
```{r}



#H2 example. 
kharkivce<-subset(ukrainepost,ukrainepost$location=="Kharkiv")

kharkiv <- get_stamenmap(bbox = c(left = 36, bottom = 49.85, 
                                  right = 36.5, top = 50.1), 
                         zoom = 13)


kharkivmap<-ggmap(kharkiv)+theme_void()
 
# Fires within bounding box
left <- 36
bottom <- 49.85
right <- 36.5
top <- 50.1

#Subset fires 
kharkivfires <- ukrfirepost[ukrfirepost$longitude >= left &
                                   ukrfirepost$longitude <= right &
                                   ukrfirepost$latitude >= bottom &
                                   ukrfirepost$latitude <= top, ]

#Event types
pr <- prop.table(table(kharkivce$event_type)) * 100
pr
```

```{r eval=FALSE, include=FALSE}
write.csv(kharkivfires,"C:/Users/Mikae/OneDrive/Dokumenter/Articles/Ukraine Submission/Appendix and Data/kharkivFIRES.csv")
write.csv(kharkivce,"C:/Users/Mikae/OneDrive/Dokumenter/Articles/Ukraine Submission/Appendix and Data/kharkivce.csv")
```

```{r}

kharkivfiresmatch <- kharkivfires[kharkivfires$acq_date %in% kharkivce$event_date, ]
kharkivce<-as.data.frame(kharkivce)



ggmap_kharkiv <- kharkivmap +
  geom_point(data = kharkivfiresmatch, aes(x = longitude, y = latitude), color = "red", size = 3) +
  geom_point(data = kharkivce, aes(x = longitude, y = latitude), color = "blue", size = 3)  +
    scalebar(x.min = 36, x.max = 36.48,
                   y.min = 49.85, y.max = 50.1,
                   dist = 5, dist_unit = "km",
                   st.bottom = FALSE, st.color = "Black",
             transform = TRUE, model = "WGS84")
 
ggmap_kharkiv 
 
```



```{r}
#Ethiopia H3
firestigray<- read.csv("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Ukraine Submission/Appendix and Data/firestigray.csv")
firesNPP<-read.csv("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Ukraine Submission/Appendix and Data/TigrayNPP.csv")
ethioregions<- read.csv("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Ukraine Submission/Appendix and Data/ethioregions.csv")
ethioacl<-read.csv("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Ukraine Submission/Appendix and Data/ethiopia.csv")

tiacl<-subset(ethioacl,ethioacl$admin1=="Tigray")

tiacl$acq_date<-as.Date(tiacl$event_date, format = "%d %B %Y")

tiacl<-subset(tiacl,tiacl$acq_date >= "2020-11-03" 
                 & tiacl$acq_date <= "2021-02-01")

unique(tiacl$event_type)

 
firestigray$acq_date<-as.Date(firestigray$acq_date, format = "%m/%d/%Y")

class(firestigray$acq_date)

firestigray<-subset(firestigray,firestigray$acq_date >= "2020-11-03" 
                 & firestigray$acq_date <= "2021-02-01")



```

```{r}
ethiopia_url <- "https://data.humdata.org/dataset/cb58fa1f-687d-4cac-81a7-655ab1efb2d0/resource/63c4a9af-53a7-455b-a4d2-adcc22b48d28/download/eth_adm_csa_bofedb_2021_shp.zip"
download.file(ethiopia_url, basename(ethiopia_url), mode = "wb")
unzip("eth_adm_csa_bofedb_2021_shp.zip")
ethiopia_poly <- readOGR("eth_admbnda_adm1_csa_bofedb_2021.shp")

tigray<-subset(ethiopia_poly,ethiopia_poly@data$ADM1_EN=="Tigray")
```

```{r}
tigray_data<-fortify(tigray)
 

tigraymap<-ggplot()+
  geom_path(data=tigray_data,aes(long,lat,group=group))+
  geom_point(data=tiacl, aes(x=longitude, y=latitude), 
             color="blue", size=1)+geom_point(data=firestigray, aes(x=longitude, y=latitude), 
             color="red", size=1)+geom_point(data=firesNPP, aes(x=longitude, y=latitude), 
             color="red", size=1)


tigraymap<-tigraymap+theme_map()

tigraymap
 
```   
```{r}
#Verified Detections
tigraysat<-read.csv("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Ukraine Submission/Appendix and Data/Tigrayfires.csv")

tighi<-subset(tigraysat,tigraysat$rating=="High")
tigmi<-subset(tigraysat,tigraysat$rating=="Mid")

unique(tigraysat$rating)

firesnearacled<-read.csv("C:/Users/Mikae/OneDrive/Dokumenter/Articles/Ukraine Submission/Appendix and Data/firesnearacled.csv")
firesnearacled<-subset(firesnearacled,firesnearacled$Rating!="Low")
```

```{r}
tigraymap2<-ggplot()+
  geom_path(data=tigray_data,aes(long,lat,group=group))+
  geom_point(data=tiacl, aes(x=longitude, y=latitude), 
             color="blue", size=1)+geom_point(data=tighi, aes(x=x, y=y), 
             color="red", size=1)+geom_point(data=tigmi, aes(x=x, y=y), 
             color="green", size=1)+geom_point(data=firesnearacled, aes(x=longitude, y=latitude), 
             color="orange", size=1)

tigraymap2<-tigraymap2+theme_map()

tigraymap2
 
```


```{r}
library(ggplot2)


# Your existing code
tigraymap3 <- ggplot() +
  geom_path(data = tigray_data, aes(long, lat, group = group)) +
  geom_point(data = tiacl, aes(x = longitude, y = latitude, color = "Tiacl"), size = 1) +
  geom_point(data = tighi, aes(x = x, y = y, color = "Tighi"), size = 1) +
  geom_point(data = tigmi, aes(x = x, y = y, color = "Tigmi"), size = 1) +
  theme_map()

# Add a color-coded legend
tigraymap3 <- tigraymap3 +
  scale_color_manual(values = c("Tiacl" = "blue", "Tighi" = "red", "Tigmi" = "green"),
                     labels = c("ACLED", "VIIRS - High", "VIIRS - Mid"))+
  labs(color = "")

tigraymap3


```


```{r}
table(tigraysat$rating)

prop.table(table(tiacl$geo_precision))


```

```{r}
write.csv(tiacl,"C:/Users/Mikae/OneDrive/Dokumenter/Articles/Ukraine Submission/Appendix and Data/tiacl.csv")

```

